Discrete model for laser driven etching and microstructuring of metallic surfaces 
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We present a unidimensional discrete solid-on-solid model evolving in time using a kinetic Monte 
Carlo method to simulate micro-structuring of kerfs on metallic surfaces by means of laser-induced 
jet-chemical etching. The precise control of the passivation layer achieved by this technique is 
responsible for the high resolution of the structures. However, within a certain range of experimental 
parameters, the microstructuring of kerfs on stainless steel surfaces with a solution of H3PO4 shows 
periodic ripples, which are considered to originate from an intrinsic dynamics. The model mimics 
a few of the various physical and chemical processes involved and within certain parameter ranges 
reproduces some morphological aspects of the structures, in particular ripple regimes. We analyze 
the range of values of laser beam power for the appearance of ripples in both experimental and 
simulated kerfs. The discrete model is an extension of one that has been used previously in the 
context of ion sputtering and is related to a noisy version of the Kuramoto-Sivashinsky equation 
used extensively in the field of pattern formation. 

PACS numbers: 61.82.Bg, 81.65.Cf, 05.10.Ln, 47.54.+r 
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I. INTRODUCTION 

Since the 1980s, laser induced wet chemical etching in 
silicon, ceramics, and metals has been an intensively used 
micro-structuring technique Q. Controlling the quality 
of the final structures is the main issue from the exper- 
imental and theoretical point of view. In a variation of 
the experimental technique developed by Metev, Stephen 
and collaborators (see Refs. and 0j)and currently im- 
plemented by Rabbow et al. (see Ref. Q), a focused laser 
beam induces an etching reaction enhanced by a coaxial 
jet of etchant producing holes and kerfs on metallic sam- 
ples within a micrometer scale. 

This technique is called laser induced jet- chemical etch- 
ing (LJE) and has been successful in fabricating supere- 
lastic microgrippers of nickel-titanium alloy. However, 
within a certain range of parameters, in the microstruc- 
turing of kerfs on stainless steel surfaces with a solution 
of H3PO4, unwanted periodic ripples appear. In this pa- 
per, it is considered that these ripples are intrinsically 
generated and belong to a wide universality class of pat- 
tern formation phenomena that emerge, for example, in 
ripple structures formed by wind over a sand bed, ion 
sputtering of various surfaces 0, @ and abrasive water- 
jet cutting Q. 

In the context of ion sputtering, Cuerno, Makse, Tom- 
massone, Harrington and Stanley (CMTHS) proposed a 
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stochastic one dimensional (ID) discrete solid-on-solid 
(SOS) model based on the competition between erosion 
and surface diffusion processes in which ripples appear at 
early stages of the evolution @ . We have extended and 
adapted the CMTHS model for the LJE case taking into 
account a moving Gaussian-distributed laser beam which 
leads to a localized heating and etching of the metallic 
sample. The motivation to use such a phenomenological 
model is originated from the limited knowledge of the in- 
teraction of diverse processes occurring in a wide range of 
scales. The dynamics of the laser light absorption, heat, 
chemical reactions, hydrodynamics, and transport phe- 
nomena is too complex to be fully modelled. While this 
extended model is quite simple, it nevertheless captures 
essential physical effects of the process, such as the unsta- 
ble temperature driven etching and the stabilizing mech- 
anism of surface diffusion. In addition, since the model 
is evolved in time using a kinetic Monte Carlo method, 
fluctuations and rough surfaces are produced naturally. 

This paper is organized as follows. In Sec. |H] we 
present the basics of the experimental setup and describe 
qualitatively some of the relevant microscopic processes 
occurring during etching. The appearance of a ripple 
regime in a series of kerfs structured with increasing laser 
power is analyzed in Sec. IIIII A hypothesis about the 
intrinsic nature of the ripples is proposed there. In Sec. 
II VI we describe and justify the extension of the CMTHS 
model analyzing the scaling properties and ripple regimes 
for uniform erosion. The application of the extended 
model to simulate the LJE is presented in Sec. |VJ Vary- 
ing the simulated laser power, a ripple regime is obtained 
and compared with the experimental one. In the conclu- 
sions we summarize achievements and shortcomings of 
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FIG. 1: (Color online) Schematic diagram of the experimental 
setup, (a) The etching cell. A focused laser induces an etching 
reaction enhanced by a coaxial running jet of etchant. (b) The 
whole experimental setup. The reaction can be observed with 
a charged-coupled device (CCD) camera and the intensity of 
the reflected light is measured with a photoresistance. Due 
to the laser induced etching the workpiece changes its state 
locally from passive to active. Changes of the potential are 
detected with an Ag/AgCl reference electrode. 



the model and suggest future improvements. 



II. LJE TECHNIQUE 
A. Experimental setup 

A schematic diagram of the etching cell is shown in 
Fig. ^a) . In this implementation of the technique, foils 
of stainless steel (Fe/Crl8/Nil0) are immersed horizon- 
tally in a solution of etchant based on 5-M H3PO4. The 
microstructuring of the samples is achieved by an argon 
ion cw-laser beam at 514nm embedded coaxially in a jet 
of liquid etchant which is directed perpendicularly to the 
surface. The intensity profile of the focused beam is as- 
sumed to be Gaussian with an estimated standard de- 
viation of a 2 //m. The laser spot diameter for the 
experiment is defined as d exp = Aa ~ 8 /im. The laser 
induced etching leads to dissolution of the metal and for- 
mation of hydrogen. Simultaneously, changes of the state 
of the electrode from passive to active produce an electro- 
chemical potential and its temporal evolution E(t) can be 
measured against an Ag/AgCl reference electrode, which 
is immersed in the etching reservoir. 

The etching cell is mounted on a computer-controlled 
mobile basis which allows us to move the sample in the 
xy-plane with respect to the laser beam with different 
feed velocities Vf [see Fig. H^b)]. The setup is automated 
allowing us to structure holes or kerfs varying the most 
relevant external parameters: laser power, etchant jet ve- 
locity, feed velocity, and etchant concentration. Details of 
the process can be observed with a CCD camera located 
above the etching cell in the same axis of the laser beam. 
Alternatively, by means of a photoresistance located in 
the same position, measurements of the intensity of the 
reflected light can be used together with the electrochem- 
ical potential for monitoring the etching dynamics. 



B. Microscopic description of the etching process 

Under normal conditions of temperature, the layer in 
contact with the liquid is passivated spontaneously thus 
isolating the metallic sample from the etchant action. 
The area below the laser spot is heated and almost imme- 
diately the heat spreads to a wider zone. Above a tem- 
perature threshold the passivation layer is removed and 
thermally activated chemical etching starts there form- 
ing the etching front which structures the kerf on the 
surface at velocity Vf. The protons of the phosphoric 
acid (H3PO4) react with the iron, nickel, and chromium 
of the steel producing hydrogen and dissolution of metal 
ions. The basic reactions can be described as 



Reduction (formation of hydrogen) : 2H + + 2e 



H 2 T 
(1) 



Oxidation (ionization of the metal) : Fe — > Fe 2+ + 2c 

with similar reactions for the ionization of the nickel and 
chromium. The etching reaction is exothermic in nature. 

When the etching reaction dissolves the metal and con- 
sumes the etchant, a thin layer of solution in contact with 
the metallic surface develops a concentration gradient 
ranging from zero on the surface to the value of the bulk 
concentration. This is called the Nernst diffusion layer 
(NDL) and within its thickness 5, the transport of ions 
of etchant and reaction products of the reactions occurs 
exclusively by diffusion, limiting the etching reaction . 
Outside this layer, convective transport maintains the 
concentration uniform at the bulk concentration. The 
value of S depends on the hydrodynamic conditions im- 
posed by the etchant jet. 

The NDL is not directly modelled in this work. We 
assume that its thickness S is almost negligible compared 
with the dimensions of the simulated topographies. It is 
worthwhile to note that the role of the Nernst diffusion 
layer is indeed essential in the real etching process. Vari- 
ations of the dynamics of the jet due to its interaction 
with the surface can affect the value of the thickness S 
and in consequence, the transport of ions and reaction 
products. For example, if in a trough of a ripple the tur- 
bulence of the etchant flow allows a growing thickness 5, 
this would imply a strong inhibition of the etching rate 
after a certain depth threshold. 

In experiments of wet-chemical etching of ceramics 
(without etchant jet), Lu et al. have found that the 
diffusive transport in the NDL is the main limiting fac- 
tor of the etching rate [lj} ■ Following their analysis and 
assuming that the etching rate decreases exponentially 
with the depth, an expression for the dependency of the 
etched depth with the feed velocity has been found to be 
in good agreement with experimental data for the LJE 
experiment using 5 — M H3PO4 (see Ref. Q). The func- 
tion of the etchant jet is to enhance the etching rates 
reducing the NDL thickness 6, to provide fresh etchant 
and to remove dissolved material and reaction products. 
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In addition, the jet creates a cooling effect that main- 
tains the heating effect of the laser spot concentrated in 
a small region. Therefore the de-passivated zone and the 
resulting etching are highly localized. When the etchant 
jet velocity is increased the etching reaction is favored 
due to more fresh etchant but on the other hand it is 
inhibited due to the cooling effect. After the laser jet 
leaves a region the surface is repassivated due to the de- 
crease of the temperature and its topography will remain 
unaltered. 



III. OBTAINED KERFS AND INTRINSIC 
RIPPLE FORMATION 

In most cases the etching front is stationary and the ob- 
tained kerf has, apart from small fluctuations, a defined 
shape with almost constant width and depth. The bot- 
tom and walls of the kerfs are rough due to the stochastic 
character of the etching reactions. On the other hand, 
for certain parameters ranges of feed velocity, laser beam 
power, or etchant jet velocity, oscillations of the etching 
front producing periodic ripples have been reported [4|. 
Here we will analyze the series of experiments with in- 
creasing power and a constant feed velocity Vf = 6 /zm/s 
shown in Fig. [3 

There are two external sources which produce periodic 
structures that can be detected in the measurements for 
stationary etching fronts. First, the peristaltic pump in- 
troduces vibrations in the etchant jet in the range of 1- 
3Hz that can be detected in the power spectra of the elec- 
trochemical potential time series, which correspond to 
periodic surface structures smaller than 10 //m. Second, 
the xy stage has an internal mechanism that controls its 
position every 40 /zm, resulting in small oscillations in the 
feed velocity Vf around a mean value. Nevertheless, the 
frequencies corresponding to the examined ripples are in- 
commensurable with the frequencies of both the etchant 
pump and the controlling mechanism of the xy stage. In 
consequence these can be discarded as external triggers 
of the obtained ripples. 

From the point of view of the theory of pattern for- 
mation in continuous media, the appearance of ripples 
originated from instabilities is not unexpected because 
the etching front is formed in a system far from equi- 
librium due to the continuous and combined action of 
the laser beam and the etchant jet. In general, one class 
of mechanisms for such instabilities arises from the exis- 
tence of constraints and conservation laws in the system 
[Tl|. In the case of water jet cutting, at high feed rates 
intrinsically generated ripples and striation patterns of 
the order of the water jet diameter are formed at the 
side walls of the cut (see Ref. [l2|). This fact has been 
used there as a criterion to distinguish between ripples 
triggered externally from the intrinsic ones. 

The diameter of the etching front is not always con- 
stant for the LJE experiment. Certainly, more delivered 
power means higher temperatures and, although the laser 
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FIG. 2: Reproduced from Rabbow y. Optical microscope 
images showing a series of kerfs structured with increasing 
powers from 250 to 650 mW (etchant 5-M H3PO4, feed veloc- 
ity 6 jim/s, etchant jet velocity 190 cm/s). The kerfs for 250, 
300, and 350 mW are the result of a uniform etching front 
which widens with power. For powers greater than 350 mW, 
an instability in the etching front comes in, and the obtained 
ripples seem to be product of a thermal runaway. A hypothe- 
sis on the intrinsic character of the ripples and some remarks 
on the mechanism that is creating this pulsating etching front 
are discussed in the text. 



spot diameter remains the same, the etching front be- 
come wider. For laser powers 250, 300, and 350 mW 
an approximately stationary etching front produces kerfs 
with increasing width and depth (see first three kerfs in 
Fig. 0). 

For powers larger than 350 mW an instability appears 
producing ripples with increasing length and width. The 
ripples look like a product of thermal runaways which 
reach much broader areas than in the stationary etching 
front case. Each thermal runaway seems to appear above 
a temperature threshold, producing a quick broadening 
of the etching front. The broadening stops when the ac- 
celerated consumption of etchant inhibits further etching 
rates and the etching front shrinks until the conditions 
for the onset of the next thermal runaway are fulfilled. 
Note that during the ripple formation the absorption of 
laser radiation is also changing according to variations 
of the slope of the etching front. Probably this has a 
reinforcing effect on the pulsating etching front. 

Based on these facts we formulate the hypothesis that 
above a power threshold an instability creates ripples. 
The resulting wavelength is in the order of magnitude of 
the diameter of the stationary etching front that would 
appear if the mechanism that causes the instability would 
not act. This could explain why the ripple length in- 
creases with the power. For powers larger than 600 mW 
the ripples seem to overlap as the kerfs become deeper. 

IV. DISCRETE MODEL 

Instead of formulating a microscopical model with the 
interaction of all possible chemical and physical pro- 
cesses, we propose a minimal phenomenological descrip- 
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(a) 



(b) 
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FIG. 3: Diagram representing the surface and the erosive and 
diffusive actions, (a), (b): sites submitted to the erosion rule. 
The probability of it being eroded is larger for the "valley" 
in (b) than for the "peak" in (a) according to an estimation 
of the curvature described in Sec. 11V A II (c), (d): diffusive 
movements. The probability of surface diffusion in (c) is close 
to f , while for (d) it is very small according to a mechanism for 
creating a positive surface tension described in Sec. IIV A 21 



tion of the etching process based on some ideas of the 
theory of far from equilibrium evolving surfaces. A num- 
ber of discrete growth models and continuum stochas- 
tic equations have been proposed to describe the kinetic 
roughening properties of surface growth and erosion |13| . 
For simulating the structures produced by the LJE tech- 
nique, we have extended and adapted a stochastic ID 
discrete solid-on-solid (SOS) model proposed by Cuerno 
et al. 8] for the evolution of ion sputtered surfaces 
(CMTHS model). Within that framework they explain 
two stages: an early time regime characterized by rip- 
ples and late time regime where the roughening shows 
self-affine scaling behavior. It has been shown that this 
model is related with a noisy version of the Kuramoto- 
Sivashinsky equation, which has been used extensively in 
the theory of unstable pattern formation |l4|. 

A remarkable feature in experiments of ion sputter- 
ing is the presence of nearly periodic ripples, aligned 
parallel or perpendicular to the bombarding ion beam 
[H, . Relating the energy of the ion beam with sput- 
tering yield, Bradley and Harper (BH) found that 
the dependency of the erosion rate with the local sur- 
face curvature induces an instability, which is responsible 
for the formation of periodic ripples with a characteris- 
tic length. Troughs are eroded faster than peaks and 
this effect can be considered as a "negative surface ten- 
sion," which competes with the smoothing mechanism of 
thermally activated surface diffusion (which is a positive 
surface tension). 

Our extended model is based on similar mechanism 
and the justification for applying it to ripples found in 
LJE is based on experimental evidence that shows that 
troughs are preferentially etched as compared to peaks. 
This can be explained as consequence of differences in 
the absorption of laser energy and the resulting heat pro- 
cesses. First, when the laser beam is acting on a trough, 
due to the geometry the reflected rays converge to the 
zone above the trough and eventually can produce mul- 



tiple reflections inside. This means that the trough and 
its neighborhood can effectively absorb more laser radi- 
ation, more heat is produced, and in consequence the 
etching rate is enhanced. In the case of peaks reflected 
rays are dispersed in all directions and there are no sec- 
ondary reflections. On the other hand, considering the 
cooling effect of the etchant jet due to heat convection, it 
is easy to imagine that peaks are cooled more efficiently 
than troughs where eddies and even stagnation of the 
fluid are more probable to appear. Then, the preferen- 
tial heating and etching within troughs result in further 
increase of the local curvature, which in turn enhances 
the secondary reflections and heating due to poorer con- 
vective heat transport by the etchant flow. 



A. Extension of the CMTHS model 

The material to be eroded is represented in 1+1D by 
a lattice composed of cells of width a and height b, and 
the surface is represented by the integer valued height hi 
where i — 1, . . . , L (see Fig. [3J . The system size L is the 
number of cells in the horizontal direction and periodic 
boundary conditions apply for the height hi. For each 
column, all the sites below the surface are occupied with 
cells, whereas all the sites above are empty. Overhangs 
are not allowed in this kind of model for simplifying both 
the simulation and the analytical approach. 

The temporal evolution of this virtual surface takes 
into account rules for representing the erosion and the 
surface diffusion processes. In this kind of kinetic Monte 
Carlo simulations each erosion or diffusion event appears 
with a rate that has to be guessed through the use of all 
the available experimental and theoretical information 
[l7|. The program defines a flat surface (one of the pos- 
sible initial conditions), selects a site, and invokes with 
probability / the erosion rule and with probability 1 — / 
the surface diffusion rule. The process of selection of 
sites or rules to be applied is performed using a random 
number generator described in Ref. [Tsj . 



1. Extended erosion rule 

The erosion probability p e for a cell at the site i is 
estimated as the product p e = p K Yi. The quantity p K 
corresponds to the probability of being eroded depend- 
ing on the curvature of the surface at the site and ac- 
counts for the unstable erosion mechanism that exists in 
the physical system. The value of p K is larger for pos- 
itive curvatures than for negative ones [see Figs. Efa) 
and[3{b)]. In the CMTHS model for ion sputtering the 
dependency of the erosion rate on the angle of incidence 
<p between the beam and a tilted portion of the surface 
is described by the sputtering yield function (from Ref. 
0): 
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ity p K ,min- The obtained values of the curvature are 
mapped by a linear transformation in such a way that for 
Ki = K max the curvature dependent erosion probability is 
Pk = Vn,max — 1 and when Ki — K m i n then p K = p K)m i n . 
In the case that the computed curvature Ki were larger 



than K max the algorithm assigns Ki 
when Ki <C K m i n then Ki — K m i n . 



Similarly 



FIG. 4: (Color online) (a) The sputtering yield dependence 
with the angle of incidence. The function used in the sim- 
ulations is Y((p) = 0.5 + 0.979V3 2 - 0.479</ and Y(0°) = 
0,y(57.3°) = 1,F(90°) = 0. (b) Absorptivity of polarized 
light versus incidence angle for a flat iron surface (based on 
Ref. llSBI. Plane of polarization parallel to the incidence 
plane (continuous line), circular polarization (dotted line), 
plane of polarization parallel to the incidence plane (dashed 
line). The coefficient of refraction is n = 3.81 and the atten- 
uation coefficient is k = 4.44. 



where ipi is the incidence angle formed by the incoming 
beam and the normal direction defined on the surface at 
the site i. We can use the same function for the LJE case 
based of the fact that absorption of polarized light by 
flat metallic surfaces has similar functional dependence 
(in the case of electric field parallel to the plane of in- 
cidence or even circular polarization, see Fig. QJ. The 
nonlinearity introduced by this yield (absorption) func- 
tion becomes relevant at late regimes when large slopes 
develop, then the ripples will be strongly distorted and 
the surface will have a rough morphology. 

We propose to replace the "box rule" described in Ref. 
@ by a direct estimation of the angles and the curva- 
tures based on a finite central differences method around 
a selected site. The first derivative or gradient of the 
surface at a point i is V» = (hi-i — hi+x)/2a and the 
angle ifi that the surface form at this site is estimated by 
ipi = arctan(Vi). This angle corresponds to the incidence 
angle used in the formula ©. The second derivative is 
V 2 = (hi-i — 2hi + hi + i)/a 2 . The curvature can be esti- 
mated using the standard formula Ki — Vf [l+(Vi) 2 ]~ 3 / 2 . 

Due to the discreteness of the height hi the values ob- 
tained with these formula vary drastically from one site to 
the other. To attenuate this problem, the values for the 
angles and curvatures are computed not only for the site 
i but also for its neighbors i — 1 and i + l. Then the mean 
value of the angle for the site i is Tpl = ((fi^i+tpi+ipi+i) /3 
and the mean curvature is: ~kI — (m—i + Ki + «i+i)/3. 
This procedure takes into account the values of five sites 
(hi-2, hi-i, hi, hi + i, hi+2) an d provides a smoothed es- 
timation of the angles and curvatures which will influence 
strongly the evolution of the topography of the surface. 

Taking into account these modifications in the algo- 
rithm, it is necessary to introduce two new parameters 
for estimating the curvature dependent erosion proba- 
bility p K . First, the maximum of the positive curva- 
ture K max (the minimum K m i n is the negative of this 
value) and second, the minimum of the erosion probabil- 



2. Surface diffusion rule 

This rule is implemented in the same way as proposed 
in the original CMTHS model. The probability of a diffu- 
sive movement of a selected cell i is evaluated selecting at 
random one of the two nearest nei ghb ors and computing 
the hopping probability (see Ref. 
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1 + exp[AHi^i±i/(k B T)Y 



(4) 



where AHi^i±i is the energy difference between the fi- 
nal and initial states of the move, Ub is the Boltzmann 
constant, and T is the temperature. This surface energy 
is defined by the Hamiltonian: 



b 2 



h-i+i) 2 , 



(5) 



where J is a coupling constant through which the nearest 
neighbor sites interact and b is the height of the cells. 
Diffusion movements that produce final states with lower 
surface energy are then highly preferred [see Figs. Etc) 
and Eld)]. 



B. Scaling of the extended model 

As usual in the field of far from equilibrium interfaces, 
we are interested in the temporal behavior of the surface 
width defined as 



W(t) = 



\ 



(6) 



where L is the system size and (h(t)) is the mean value 
of the all heights of the surface at the time t. The 
width W(t) can be considered as a characterization of the 
roughness and it is evaluated averaging over a number of 
different realizations of the random number seed. In or- 
der to evaluate an eventual scaling behavior W(t) ~ t 13 
within an interval, the growth exponent [3 is estimated 
by the method of the consecutive slopes (see Ref. p~3| . 
p. 305). For all the figures and analyses of this section, 
the time unit t is chosen to correspond to L erosion or 
diffusion rule invocations. 

In order to verify the connection of the extended model 
with the original CMTHS model, the evolution of the 
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FIG. 5: (Color online) Comparison of the scaling of the sur- 
face width for the CMTHS and extended models. For the 
CMTHS model: The surface width W c (t) (bold continuous 
line, scale on the left) and its growth exponent /3 C (bold dashed 
line, scale on the right). Parameters: L — 2048, / = 0.5, 
J/(k B T) = 5. For the extended model: The width W c (t) 
(thin continuous line, scale on the left) and its growth ex- 
ponent f3 c (thin dashed line, scale on the right). The extra 
parameters of the extended model are the closest possible to 
the CMTHS model: a = 1, b = 1, p K:min = 1/7, K max = 2. 



surface width for both is compared in Fig. [SJ The system 
size is L = 2048 and the width values were averaged over 
100 different realizations. The probability of invoking 
the erosion rule is / = 0.5 (the diffusion rule is invoked 
with probability 1 — / = 0.5) and the constant associated 
with the diffusion is J /{ksT) — 5. The extra parameters 
of the extended model (a = 1, b = 1, p Kl min = 1/7, 
Kmax = 2) have been chosen to reproduce as closely as 
possible the box rule for the erosion probability of the 
CMTHS model. For both cases the same "yield" function 
is used Y((p) = 0.5 + 0.97V - 0.479(^ 4 . 

The scaling properties are similar, showing first a 
rough interface at early times, then a strong increase of 
the growth exponent (3 due to the onset of the instabil- 
ity that creates ripples, and finally a drop due to the 
stabilizing and roughening effect of the non-linear terms. 
However, the scaling behavior is not identical. A precise 
identification of the limits of the ripple regimes is dif- 
ficult because they depend on each realization and the 
criterion to distinguish between fluctuations and proper 
ripples. The ripple regime for the CMTHS model can be 
estimated to be t ss (300,1000); while for the extended 
model, the onset of the instability occurs earlier and the 
ripple regime lasts longer. The ripples in the extended 
model are larger in amplitude and present a quasi sinu- 
soidal shape. 

The application of the extended model for the sim- 
ulation of the LJE that is presented in the next section 
requires to exploit the flexibility of the new parameters to 
generate ripples with a characteristic length in the order 
of the size of active etching front. In addition, the ripples 
should be clearly distinguishable from the inherent fluc- 
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FIG. 6: (Color online) Temporal evolution of the interface 
width W c (t) (continuous line, scale on the left) and growth 
exponent j3 e (dashed line, scale on the right) of the extended 
model for the parameter set used in Sec. [V] The surface's 
morphology at times corresponding to t = 51 (o symbol), 
t = 442 (□ symbol), t = 11437 (A symbol) and t = 
10 (0 symbol) is analyzed in Fig. |7| Extended model param- 
eters: Y(ip) = 1, L = 2048, / = 0.045, J/(k B T) = 1, a = 20, 

6=1, PK.min = 0.1, K max = 0.0004. 



tuations and roughness produced by the stochastic char- 
acter of the model. We will apply the erosion and surface 
diffusion rules within a moving probability distribution 
that is proportional to a temperature field generated by 
the absorption of light on the region illuminated by the 
focused laser. Details will be explained in Sec. The 
active etching zone spans a region wider than the laser 
spot. We will assume that variations introduced by the 
angle dependence of the absorption of light at the laser 
spot expressed by the yield (absorption) function Y(ip) 
can be neglected in a first approximation. 

Accordingly, in the rest of this section, we will an- 
alyze the scaling properties and the morphology evo- 
lution using Y(tp) = 1 and the same parameters that 
will be used for the simulations of the LJE presented 
in Sec. The probability of invoking the erosion rule 
is / = 0.045 (the diffusion rule is invoked with proba- 
bility 1 — / = 0.955) and the constant associated with 
the diffusion is J/{ksT) = 1. For the curvature depen- 
dent erosion probability p K the values a = 20, 6=1, 
PK.min = 0.1, K m ax — 0.0004 have been used. The sys- 
tem size is L — 2048 and the width values were averaged 
over 30 different realizations. 

Figure [5] shows the temporal evolution of the width 
W e (t) and its growth exponent (3 e for this parameter set. 
The ripple regime persists much longer than in the cases 
presented in Fig. [SJand the obtained ripples have a larger 
amplitude. For times t < 150 the exponent is lower 
than the value 0.5 characteristic for random erosion. The 
strong increase of the value of f3 after t > 100 is associ- 
ated with the onset of the ripple regime. During the 
interval 10 3 < t < 10 4 there is an oscillation of (3 but the 
values remain relatively high. The origin and meaning of 
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FIG. 7: (Color online) Four different stages of the evolution 
of the surface for the extended model corresponding to the 
times referred to in Fig. |S| The left column shows the height 
h* (note that all vertical scales are different). The real system 
size is L = 2048 but only 300 points are shown in order to 
appreciate details of the morphology of the surface. The right 
column shows log-log plots of the corresponding ensemble av- 
erage of the power spectral density EA-PSD (lower curve) and 
the smoothed ensemble average of the power spectral density 
SEA-PSD (upper curve, shifted in the vertical direction for 
clarity of presentation). The local maxima of the SEA-PSD 
[indicated by the (V) symbols] is used to estimate the mean 
value of the ripple length / of each stage. 



this oscillation is currently unknown. After t w 10 4 the 
growth exponent (3 shows a slow increase. 

Figure shows a portion of the surface for the four 
times that are indicated with symbols in Fig. [fJl The left 
column shows the height h*, which is the height minus 
its average value at each time h* = hi — (hi). The right 
column shows the corresponding ensemble average of the 
power spectral density EA-PSD computed with 30 real- 
izations and L = 2048. The PSD is represented in a loga- 
rithmic scale with arbitrary units and the horizontal axis 
corresponds to the wave number k. A moving window 
average over ten points is applied on the spectra and the 
resulting smoothed ensemble average of the power spec- 
tral density SEA-PSD appears above the EA-PSD curve. 
The local maxima of the SEA-PSD (indicated by the V 
symbols) can be used to estimate the mean value of the 
ripple length I of each stage. 

In the rough surface corresponding to t = 51 
(o symbol) the integer values of the heights are visible. 
For t — 442 (□ symbol) the ripples start to appear but 
their shape and length are highly irregular. More soft and 
larger ripples are found for times in the regime around 



to t = 11437 (A symbol) and the length estimated by 
the SEA-PSD method is approximately I = 23a. It is ob- 
served that the ripple length is increasing with time due 
to the merging of ripples: small ripples are assimilated 
by contiguous larger ripples, which in turn develop a si- 
nusoidal shape. For time t = 10 5 (0 symbol) the ripple 
length is approximately 40a. 

The increasing ripple length is a deviation from lin- 
ear analysis predictions on the Kuramoto-Sivashinsky de- 
scription. It is originated from nonlinear terms that can 
be examined in the derivation of the continuum equa- 
tion for this discrete model as it is shown in Ref. |21| . 
Analogous coarsening phenomena have been reported for 
experiments on ion sputtering [23, H3 and laser ablation 
[24]. among others. 

V. SIMULATION OF KERFS FORMATION IN 

LJE 

In order to simulate the joint action of the laser 
beam and the etchant jet it is necessary to estimate the 
heat spreading resulting from the absorption of the laser 
beam. The intensity of the moving laser can be expressed 
by means of 

G( r ,t) = — S=e-< r -" ^l^ 2 \ (7) 
cry 27r 

where r is the spatial coordinate, t the time, a corre- 
sponds to the standard deviation, and Vf is the feed 
velocity. Given the thickness of the metallic sample 
(~ 200 fim) and the relatively large dwell times of the 
laser beam, the back surface of the sample reaches ap- 
proximately the same temperature as the front surface 
on which the radiation is incident. For this "thermally 
thin approximation" (see Chapter 3 in Ref. 25]), the 
computation of the temperature field requires numerical 
integration that will be reported in detail elsewhere. Re- 
gardless of variations related to the choice of the various 
thermal constants and parameters, a generic profile of 
the temperature field can be identified. The three curves 
shown in Fig. [SJare proportional to such a profile that can 
be described to be similar to a Gaussian shape within the 
laser spot while for larger distances it decays as —log{r). 

We define the etching probability distribution Pe of 
applying the erosion and surface diffusion rules to be di- 
rectly proportional to the estimated temperature field. 
Introducing an additional factor between and 1 to the 
Pe distribution, it is possible to simulate laser beam pow- 
ers between and 100 %. In order to simulate the fact 
that etching occurs only above a certain temperature, a 
depassivation threshold is introduced and Pe is defined 
to be zero below it (see Fig. |HJ. Therefore the etching 
probability distribution has a finite diameter d, which 
depends on power and depassivation threshold. Figure |S| 
shows Pe for three different power levels 20, 60 and 100 % 
and depassivation threshold of 0.09. The diameter and 
amplitude increase proportionally to the power. In order 
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FIG. 8: (Color online) Etching probability distribution Pe 
associated with the estimation of the temperature field pro- 
duced by a Gaussian beam with standard deviation a — 2a 
for beam powers: 20, 60, and 100 %. The etching probability 
Pe determines how frequently the erosion and surface diffu- 
sion rules are applied, and is used to simulate the joint action 
of the moving laser beam and etching jet for different powers. 
The Gaussian distribution depicted in the figure has the same 
standard deviation of the laser beam. The vertical lines are 
located at r = ±2a. 



to compare the diameters d of the three Pe distributions 
(which are 16a, 36a, and 46a, respectively) with the laser 
spot diameter used in the simulations (d SJ ; m = 4er = 8a), 
the corresponding Gaussian distribution is also shown in 
the figure. 

In the simulation the surface is kept fixed while the 
laser beam and the etching probability distribution move 
with velocity v / . Within the etching front a combined ac- 
tion of erosion and surface diffusion rules occurs during a 
characteristic dwell time d/vf. Depending upon Vf and 
the power, the etching front will eventually develop insta- 
bilities giving rise to ripple structures. After the etching 
front leaves the region the surface has acquired a topog- 
raphy that, depending on model parameters, corresponds 
roughly to one of the stages analyzed in Fig. El 

Figure shows the creation process of a single ripple 
for a feed velocity 2 x 10 _6 a/step and etching probability 
distribution Pe corresponding to power 60 % shown in 
Fig. 03 The extended model parameters are the same as 
used in Figs. Eland0of the previous section. Typically, 
a valley is created at the forefront of the moving Pe from 
a small but growing depression on the surface. When the 
center of Pe passes through a valley, the rate of erosion 
increases, and the valley grows as long as the rearmost 
part of the Pe is acting. The peaks between the valleys 
are eroded at lower rate due to their negative curvature. 
In summary, the local and temporal action of the etching 
probability distribution works as an amplifier of small 
instabilities on the surface. 

Figure^Jshows the profiles of the kerfs structured with 
etching probability distributions corresponding to pow- 
ers ranging from 20 to 100 % and constant feed velocity. 
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FIG. 9: (Color online) Four stages of the formation of a single 
ripple. The standard deviation of the moving Gaussian beam 
is a — 2a. Its center is represented by the vertical dashed 
line. The vertical continuous lines are located at ±2cr in a 
comoving system. The etching probability distribution Pe 
corresponding to 60 % of power and depassivation threshold 
0.09 shown in Fig. |H] is represented by the gray-scale bar at 
the bottom of each frame. A marker illustrates the evolution 
of a point on the surface at a fixed position i = 200. The 
feed velocity is 2 x 10 _6 a/step. Extended model parameters: 
Y(y) = 1, / = 0.045, J/(k B T) = 1, a = 20, 6 = 1, p K . min = 
0.1, n max = 0.0004. 



For increasing power, the diameter d of the corresponding 
Pe distributions grows leading to an increase of the dwell 
time d/vf. Together with the increase of the overall am- 
plitude of the probability Pe, the resulting topographies 
correspond to later stages of the evolution shown in Fig. 
For power 20 % the kerf is shallow and rough because 
its topography corresponds to early times of the evolution 
where the ripple regime is not yet reached. For powers 
ranging from 30 to 90 % regular ripples with increasing 
wavelength are obtained. For beam power 100%, after a 
few initial transient oscillations, the bottom of the kerf 
becomes approximately fiat with small roughness. 

The obtained ripples are very regular and periodic. An 
ensemble average over ten realizations with 16 384 points 
of the bottom of the kerf allow us to identify in the power 
spectral density a definite ripple length with negligible 
uncertainty. Figure lllf a) shows that the ripple length 
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FIG. 10: (Color online) Profiles of the kerfs for beam pow- 
ers 20, 30, . . . , 100 % structured with the temperature related 
etching probability distributions Pe shown in Fig. |H1 All kerfs 
are structured with the same feed velocity 2 x 10 -6 a/step. 
For power 20 % the kerf is shallow and rough while for pow- 
ers from 30 to 90 % a ripple regime with increasing ripple 
length appears. For beam power 100 % the bottom of the kerf 
again becomes approximately flat with small roughness. The 
extended model parameters are the same as in the previous 
figure. 



grows almost linearly with the power (o symbols). The 
diameter d of the etching probability distributions (de- 
picted with the A symbols) also grows with the power. 
In order to compare with the ripple lengths, the diam- 
eter of the laser spot d S i m = 8a used in the simulation 
is indicated in the figure with the horizontal line. For 
powers < 20 % and for powers > 95 % the obtained kerfs 
are uniform. 

In order to compare our model with the experiment, 
Fig- IllF bl shows the ripple lengths corresponding to 
powers 400, 450, 500, and 550 mW of Fig. El The rip- 
ple length increases proportionally with the laser power. 
The diameter of the laser spot d exp w 8 /im used in the 
experiment is indicated in the figure with the horizon- 
tal line. Of course, four points are not enough to draw 
conclusions about this dependency. We will explore fur- 
ther this correspondence when more experimental data 
become available. Summarizing, the simulations show a 
qualitative agreement of our hypotheses with the exper- 
iment: the etching front covers a region wider than the 
laser spot and the resulting characteristic ripple length is 
proportional to the diameter of the active etching zone. 
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FIG. 11: (Color online) (a) Simulation. Increasing ripple 
length (o symbols) for beam powers 25, 30, 35, . . . , 90 % cor- 
responding to the same parameter set as in the two previous 
figures. For powers 15 and 20 % the kerfs are shallow and 
uniform. For 95 and 100 % the kerfs become also uniform af- 
ter a few transient oscillations. The diameters d (A symbols) 
of the corresponding etching probability distributions Pe in- 
crease with the power, (b) Experiment. Ripple length (o 
symbols) of the kerfs shown in Fig. [5] corresponding to pow- 
ers 400 — 550 mW. The kerfs corresponding to powers < 350 
mW and > 600 mW are considered uniform. For both (a) 
and (b), the region between the two dotted vertical lines de- 
fine approximately the power interval where periodic ripples 
appear. 



VI. CONCLUSIONS 

The discrete extended model applied to the simula- 
tion of formation of kerfs reproduces qualitatively some 
of the main characteristics of microstructures produced 
by the LJE experiment. The basic mechanism based on 
the competition of curvature dependent erosion and sur- 
face diffusion generates a ripple regime, which appears 
between roughening regimes at earlier and later stages 
of the surface evolution. Compared with the original 
CMTHS model, the extended model produces longer time 
intervals where ripples occur, which is more appropriate 
for the simulation of the moving etching source of the LJE 
experiment. This is accomplished by applying the erosion 
and diffusion rules within a comoving etching probability 
distribution that is proportional to an estimated temper- 
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ature field produced by the laser spot. Taking advantage 
of the probabilistic nature of the Monte Carlo method, it 
is possible to simulate the laser power distribution, laser 
spot diameter and feed velocity. In analogy with the ex- 
periments, variations of the laser power reveal a regime 
with an unstable etching front and in consequence, reg- 
ular ripples in the bottom of the kerfs. In addition, the 
ripple length is growing with increasing power. Outside 
of this interval of power values, the kerfs are uniform with 
small roughness at the bottom. 

This simple phenomenological model does not pretend 
to simulate all details of the complex behavior of the pro- 
cess. Instead, it can be considered as a guide to propose 
experiments that reveal more relevant features of the rip- 
ple formation and the structuring process in general. The 
next step is to incorporate the slope dependency of the 
absorption of the laser polarized light into the model. In 
addition, including the dynamics of etchant concentra- 
tion in the Nernst diffusion layer together with an addi- 



tional heat source due to the exothermic reactions could 
allow us to understand the conditions for the onset of the 
thermal runaway responsible of ripple formation. For a 
better comparison with the experimental results a fur- 
ther obvious step is the generalization of the extended 
model to 2+1 dimensions. 
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